#Import des librairies utilisées
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(ggplot2)
library(tidyr)
library(tseries)
## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo
library(forecast)
library(TTR)
#Définition du répertoire de travail
path = "/home/utilisateur/Documents/Computer Science/ML/M1/time_series/time-series_projet/Dataset"
#Chargement du fichier day.csv et transformation en dataframe
data_day <- read.csv(file = paste(path, "/day.csv", sep=""))
df_data_day <- data.frame(data_day)
#Conversion de la variable date en format date
df_data_day$dteday <- as.Date(data_day$dteday)
#Renommage des saisons
df_data_day$season <- factor(format(df_data_day$season, format="%A"),
levels = c("1", "2","3","4"),
labels = c("Printemps", "Eté", "Automne", "Hiver"))
#Création de nouvelles variables numériques dé-normalisées pour faciliter l'interprétation
denormalize <- function(x, minval, maxval) {
return (x*(maxval-minval) + minval)}
df_data_day$true_temp <- Map(denormalize, df_data_day$temp, -8, 39)
df_data_day$true_temp <- df_data_day$temp*41
df_data_day$true_atemp <- df_data_day$atemp*50
df_data_day$true_windspeed <- df_data_day$windspeed*67
df_data_day$true_hum <- df_data_day$hum*100
#Visualisation et résumé statistique
head(df_data_day)
## instant dteday season yr mnth holiday weekday workingday weathersit
## 1 1 2011-01-01 Printemps 0 1 0 6 0 2
## 2 2 2011-01-02 Printemps 0 1 0 0 0 2
## 3 3 2011-01-03 Printemps 0 1 0 1 1 1
## 4 4 2011-01-04 Printemps 0 1 0 2 1 1
## 5 5 2011-01-05 Printemps 0 1 0 3 1 1
## 6 6 2011-01-06 Printemps 0 1 0 4 1 1
## temp atemp hum windspeed casual registered cnt true_temp
## 1 0.344167 0.363625 0.805833 0.1604460 331 654 985 14.110847
## 2 0.363478 0.353739 0.696087 0.2485390 131 670 801 14.902598
## 3 0.196364 0.189405 0.437273 0.2483090 120 1229 1349 8.050924
## 4 0.200000 0.212122 0.590435 0.1602960 108 1454 1562 8.200000
## 5 0.226957 0.229270 0.436957 0.1869000 82 1518 1600 9.305237
## 6 0.204348 0.233209 0.518261 0.0895652 88 1518 1606 8.378268
## true_atemp true_windspeed true_hum
## 1 18.18125 10.749882 80.5833
## 2 17.68695 16.652113 69.6087
## 3 9.47025 16.636703 43.7273
## 4 10.60610 10.739832 59.0435
## 5 11.46350 12.522300 43.6957
## 6 11.66045 6.000868 51.8261
summary(df_data_day)
## instant dteday season yr
## Min. : 1.0 Min. :2011-01-01 Printemps:181 Min. :0.0000
## 1st Qu.:183.5 1st Qu.:2011-07-02 Eté :184 1st Qu.:0.0000
## Median :366.0 Median :2012-01-01 Automne :188 Median :1.0000
## Mean :366.0 Mean :2012-01-01 Hiver :178 Mean :0.5007
## 3rd Qu.:548.5 3rd Qu.:2012-07-01 3rd Qu.:1.0000
## Max. :731.0 Max. :2012-12-31 Max. :1.0000
## mnth holiday weekday workingday
## Min. : 1.00 Min. :0.00000 Min. :0.000 Min. :0.000
## 1st Qu.: 4.00 1st Qu.:0.00000 1st Qu.:1.000 1st Qu.:0.000
## Median : 7.00 Median :0.00000 Median :3.000 Median :1.000
## Mean : 6.52 Mean :0.02873 Mean :2.997 Mean :0.684
## 3rd Qu.:10.00 3rd Qu.:0.00000 3rd Qu.:5.000 3rd Qu.:1.000
## Max. :12.00 Max. :1.00000 Max. :6.000 Max. :1.000
## weathersit temp atemp hum
## Min. :1.000 Min. :0.05913 Min. :0.07907 Min. :0.0000
## 1st Qu.:1.000 1st Qu.:0.33708 1st Qu.:0.33784 1st Qu.:0.5200
## Median :1.000 Median :0.49833 Median :0.48673 Median :0.6267
## Mean :1.395 Mean :0.49538 Mean :0.47435 Mean :0.6279
## 3rd Qu.:2.000 3rd Qu.:0.65542 3rd Qu.:0.60860 3rd Qu.:0.7302
## Max. :3.000 Max. :0.86167 Max. :0.84090 Max. :0.9725
## windspeed casual registered cnt
## Min. :0.02239 Min. : 2.0 Min. : 20 Min. : 22
## 1st Qu.:0.13495 1st Qu.: 315.5 1st Qu.:2497 1st Qu.:3152
## Median :0.18097 Median : 713.0 Median :3662 Median :4548
## Mean :0.19049 Mean : 848.2 Mean :3656 Mean :4504
## 3rd Qu.:0.23321 3rd Qu.:1096.0 3rd Qu.:4776 3rd Qu.:5956
## Max. :0.50746 Max. :3410.0 Max. :6946 Max. :8714
## true_temp true_atemp true_windspeed true_hum
## Min. : 2.424 Min. : 3.953 Min. : 1.500 Min. : 0.00
## 1st Qu.:13.820 1st Qu.:16.892 1st Qu.: 9.042 1st Qu.:52.00
## Median :20.432 Median :24.337 Median :12.125 Median :62.67
## Mean :20.311 Mean :23.718 Mean :12.763 Mean :62.79
## 3rd Qu.:26.872 3rd Qu.:30.430 3rd Qu.:15.625 3rd Qu.:73.02
## Max. :35.328 Max. :42.045 Max. :34.000 Max. :97.25
g1 <- ggplot(df_data_day, aes(x=season, y=true_temp)) +
geom_boxplot()
g1 + ggtitle("Boxplots des températures par saison") +
xlab("Saison") + ylab("Température (°C)")
La période la plus chaude correspond à l’automne et la plus froide au printemps.
mean(df_data_day$true_temp)
## [1] 20.31078
median(df_data_day$true_temp)
## [1] 20.43165
La température moyenne est de 20,3°C (0,495° normalisé) et la température médiane est de 20,4°C (0,498° normalisé).
#Visualisation de la variable count par rapport à la variable temp
g2 <- ggplot(df_data_day, aes(true_temp, cnt)) +
geom_point() +
geom_smooth()
g2 + ggtitle("Locations en fonction de la température") +
xlab("Température (°C)") + ylab("Nombre de locations")
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
#Visualisation de la variable count par rapport à la variable atemp
g3 <- ggplot(df_data_day, aes(true_atemp, cnt)) +
geom_point() +
geom_smooth()
g3 + ggtitle("Locations en fonction de la température ressentie") +
xlab("Température ressentie (°C)") + ylab("Nombre de locations")
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
#Création de la variable mean_true_temp_atemp (moyenne quotidienne de la température et de la température ressentie)
df_data_day$mean_true_temp_atemp <- rowMeans(select(df_data_day, c(true_temp, true_atemp)))
#Visualisation de la variable count par rapport à la variable mean_true_temp_atemp
g4 <- ggplot(df_data_day, aes(mean_true_temp_atemp, cnt)) +
geom_point() +
geom_smooth()
g4 + ggtitle("Locations en fonction de la température moyenne \n (réelle et ressentie)") +
xlab("Température moyenne (°C)") + ylab("Nombre de locations")
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
D’après les nuages de points il semble exister une relation positive entre le nombre de locations et la température jusqu’à ce que les températures soient trop chaudes, au delà de 30°C environ. La relation devient alors négative.
#Calcul des coefficients de corrélation
print(c("R2 locations/température:", cor(df_data_day$true_temp, df_data_day$cnt)))
## [1] "R2 locations/température:" "0.627494009033492"
print(c("R2 locations/temp. ressentie:", cor(df_data_day$atemp, df_data_day$cnt)))
## [1] "R2 locations/temp. ressentie:" "0.631065699849181"
print(c("R2 locations/temp. moyenne:", cor(df_data_day$mean_true_temp_atemp, df_data_day$cnt)))
## [1] "R2 locations/temp. moyenne:" "0.630660733761854"
Les corrélations sont toutes d’environ 0,63. Ces variables sont donc corrélées.
df_data_day %>%
group_by(data_day$mnth) %>%
summarise_at(vars(true_temp, true_hum, true_windspeed, cnt),
list(mean = mean))
## # A tibble: 12 × 5
## `data_day$mnth` true_temp_mean true_hum_mean true_windspeed_mean cnt_mean
## <int> <dbl> <dbl> <dbl> <dbl>
## 1 1 9.69 58.6 13.8 2176.
## 2 2 12.3 56.7 14.5 2655.
## 3 3 16.0 58.8 14.9 3692.
## 4 4 19.3 58.8 15.7 4485.
## 5 5 24.4 68.9 12.3 5350.
## 6 6 28.0 57.6 12.4 5772.
## 7 7 31.0 59.8 11.1 5564.
## 8 8 29.1 63.8 11.6 5664.
## 9 9 25.3 71.5 11.1 5767.
## 10 10 19.9 69.4 11.7 5199.
## 11 11 15.1 62.5 12.3 4247.
## 12 12 13.3 66.6 11.8 3404.
model_registered <- lm(registered~true_temp, data = df_data_day)
print(c("R2 for bike rentals registered:", summary(model_registered)$r.squared))
## [1] "R2 for bike rentals registered:" "0.291612923597918"
model_casual <- lm(casual~true_temp, data = df_data_day)
print(c("R2 for bike rentals casual:", summary(model_casual)$r.squared))
## [1] "R2 for bike rentals casual:" "0.295158223619129"
model_ttl <- lm(cnt~true_temp, data = df_data_day)
print(c("R2 for bike rentals total:", summary(model_ttl)$r.squared))
## [1] "R2 for bike rentals total:" "0.393748731372924"
#Affichage des résultats du modèle complet (registered et casual)
summary(model_ttl)
##
## Call:
## lm(formula = cnt ~ true_temp, data = df_data_day)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4615.3 -1134.9 -104.4 1044.3 3737.8
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1214.642 161.164 7.537 1.43e-13 ***
## true_temp 161.969 7.444 21.759 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1509 on 729 degrees of freedom
## Multiple R-squared: 0.3937, Adjusted R-squared: 0.3929
## F-statistic: 473.5 on 1 and 729 DF, p-value: < 2.2e-16
plot(model_ttl, col = "green")
Il y a un lien car la p-value est inférieure à 5% et donc l’hypothèse nulle doit être rejeté mais ce lien est faible. En effet la variable température explique moins de 40% (R2 = 0,39) des locations. D’autres variables doivent être prises en compte dans le modèle.
g5 <- ggplot(df_data_day, aes(dteday, cnt)) +
geom_point() +
geom_smooth(method=lm)
g5 + ggtitle("Locations en fonction de la date") +
xlab("Date") + ylab("Nombre de locations")
## `geom_smooth()` using formula = 'y ~ x'
Il y a une tendance générale avec une croissance entre début 2011 et fin 2012. Il y a également une saisonalité : on remarque une hausse des locations durant l’été et l’automne et une baisse durant l’hiver et le printemps.
#Recherche des outliers dans la variable cnt
boxplot.stats(df_data_day$cnt)$out
## integer(0)
#Recherche des valeurs NA et Null
sum(is.na(df_data_day))
## [1] 0
is.null(df_data_day)
## [1] FALSE
Il n’y a ni outliers, ni valeurs NA, ni valeurs Null
#Création d'un time series
ts_cnt <- ts(df_data_day$cnt, start=c(2011,1), end=c(2013,1), frequency=365)
Nous travaillons avec données quotidiennes entre le 01/01/2011 (start=c(2011,1)) et le 31/12/2012 (end=c(2013,1)), le paramètre frequency doit donc être défini sur 365 (365 jours de l’année).
#Choix de la méthode de lissage
ts_cnt_HW <- HoltWinters(ts_cnt)
La moyenne mobile ne prend pas en compte les tendances et les saisonnalités. Or nos données en comporte. La méthode de lissage exponentiel de Holt-Winters a des paramètres permettant de les prendre en compte et est donc plus adaptée.
plot(ts_cnt_HW)
plot(fitted(ts_cnt_HW))
La décomposition de la série temporelle permet d’estimer les effets de la tendance et de la saisonnalité. Sur le graphique on voit que la tendance a été retirée et la saisonalité a été atténuée.
#Affichage de l'ACF et de la PACF
ts_cnt_HW_fitted <- as.data.frame(fitted(ts_cnt_HW))
ts_cnt_HW_val <- ts_cnt_HW_fitted$xhat
ggtsdisplay(ts_cnt_HW_val)
La PACF décroit rapidement ce qui suggère un modèle MA(q). Mais sur l’ACF de trop nombreuses valeurs se suivent en dehors de l’intervalle de confiance. Il faut différencier la série car elle n’est pas stationaire.
#Vérification du nombre de différenciation nécessaire
ndiffs(ts_cnt_HW_val)
## [1] 1
Il faut différencier la série une fois.
#Différentiation
ts_cnt_HW_diff <- diff(ts_cnt_HW_val)
plot(ts_cnt_HW_diff)
Le bruit est réparti aléatoirement autour de 0.
#Nouvel affichage de l'ACF et de la PACF
ggtsdisplay(ts_cnt_HW_diff)
Cette fois l’ACF décroit rapidement, ce qui suggère un modèle AR(p) avec p = 5.
#Réalisation du test Augmented Dickey-Fuller afin d'obtenir des informations plus précises
adf.test(ts_cnt_HW_diff, alternative="stationary")
## Warning in adf.test(ts_cnt_HW_diff, alternative = "stationary"): p-value smaller
## than printed p-value
##
## Augmented Dickey-Fuller Test
##
## data: ts_cnt_HW_diff
## Dickey-Fuller = -10.459, Lag order = 7, p-value = 0.01
## alternative hypothesis: stationary
La p-value est inférieure à 0,05 ce qui indique que la série est stationaire. Après différenciation nous avons donc un modèle AR(5).
ts_cnt_HW_arima <- arima(ts_cnt_HW_diff, order=c(5,0,0))
ts_cnt_HW_arima
##
## Call:
## arima(x = ts_cnt_HW_diff, order = c(5, 0, 0))
##
## Coefficients:
## ar1 ar2 ar3 ar4 ar5 intercept
## -0.5264 -0.4173 -0.3613 -0.3099 -0.2114 3.0679
## s.e. 0.0511 0.0558 0.0569 0.0558 0.0512 18.7755
##
## sigma^2 estimated as 1018331: log likelihood = -3042.92, aic = 6099.85
Avec 5 paramètres, ce modèle est complexe. Il faut en chercher un autre qui l’est moins.
#Retrait de l'effet de la saison
ts_cnt_stl <- stl(ts_cnt,"periodic")
ts_cnt_season_adj <- seasadj(ts_cnt_stl)
plot(ts_cnt_season_adj)
On constate qu’il n’y a plus de saisonnalité mais qu’il reste une tendance.
#Vérification du nombre de différenciation nécessaire
ndiffs(ts_cnt_season_adj)
## [1] 1
Il faut différencier la série une fois.
#Différentiation
ts_cnt_season_adj_diff <- diff(ts_cnt_season_adj)
plot(ts_cnt_season_adj_diff)
L’effet de la tendance a été retiré.
#Nouveau test Augmented Dickey-Fuller
adf.test(ts_cnt_season_adj_diff)
## Warning in adf.test(ts_cnt_season_adj_diff): p-value smaller than printed p-
## value
##
## Augmented Dickey-Fuller Test
##
## data: ts_cnt_season_adj_diff
## Dickey-Fuller = -13.998, Lag order = 8, p-value = 0.01
## alternative hypothesis: stationary
La p-value est plus petite que 0,05, la série est donc stationnaire.
#Auto-ARIMA
ts_cnt_auto_arima <- auto.arima(ts_cnt_season_adj)
ts_cnt_auto_arima
## Series: ts_cnt_season_adj
## ARIMA(1,1,1)
##
## Coefficients:
## ar1 ma1
## 0.2553 -0.9114
## s.e. 0.0418 0.0179
##
## sigma^2 = 401844: log likelihood = -5745.37
## AIC=11496.75 AICc=11496.78 BIC=11510.53
L’auto-arima a trouvé un modèle de paramètres 1,1,1
#Affichage des résidus
tsdisplay(ts_cnt_auto_arima$residuals, lag=20)
Il y a plusieurs valeurs hors de l’intervalle de confiance, donc qui ne sont pas du bruit blanc et non expliquée par le modèle.
#Test de Shapiro
shapiro.test(ts_cnt_auto_arima$residuals)
##
## Shapiro-Wilk normality test
##
## data: ts_cnt_auto_arima$residuals
## W = 0.95282, p-value = 1.515e-14
En effet, le résultat du test de Shapiro indique que les résidus ne suivent pas une distribution normale (la p value est plus petite que 0,05). Il reste donc de l’information non traitée par le modèle.
#Test de Ljung-Box
Box.test(ts_cnt_auto_arima$residuals, type="Ljung-Box")
##
## Box-Ljung test
##
## data: ts_cnt_auto_arima$residuals
## X-squared = 0.0070763, df = 1, p-value = 0.933
Toutefois d’après le test de Ljung-Box, les résidus ne sont pas auto-corrélés (la p-value est plus grande que 0,05).
Le modèle auto-ARIMA ayant peu de paramètres, un modèle plus complexe pourrait expliquer d’avantage les informations restées dans le bruit. Nous allons évaluer 100 modèles en testant différentes valeurs pour les paramètres p et q afin de rechercher le modèle qui minimise le score AIC.
#Itération sur les paramètres p et q entre 0 et 9
aic.values <- c()
for (p in (0:9)){
for (q in (0:9)){
ts_cnt_season_adj_arima <- arima(ts_cnt_season_adj, order = c(p,1,q), method="ML")
aic.values <- c(aic.values, ts_cnt_season_adj_arima$aic)
}
}
#Modèle ayant le plus petit AIC
print(c("Modèle avec le plus petit AIC:", which.min(aic.values)))
## [1] "Modèle avec le plus petit AIC:" "98"
#Valeur de son AIC
print(c("Valeur de l'AIC la plus basse:", aic.values[98]))
## [1] "Valeur de l'AIC la plus basse:" "11431.6205683934"
Le modèle qui minimise l’AIC est donc un modèle ARIMA d’ordre (9,1,7).
#Création du modèle de paramètres p = 9 et q = 8.
ts_cnt_season_adj_arima <- arima(ts_cnt_season_adj, order = c(9,1,7), method="ML")
## Warning in log(s2): NaNs produced
## Warning in log(s2): NaNs produced
## Warning in log(s2): NaNs produced
## Warning in arima(ts_cnt_season_adj, order = c(9, 1, 7), method = "ML"): possible
## convergence problem: optim gave code = 1
#Affichage du modèle
ts_cnt_season_adj_arima
##
## Call:
## arima(x = ts_cnt_season_adj, order = c(9, 1, 7), method = "ML")
##
## Coefficients:
## ar1 ar2 ar3 ar4 ar5 ar6 ar7 ar8
## -0.3747 0.1433 -1.1393 -0.446 -0.0106 -0.1185 -0.0272 0.0441
## s.e. 0.1686 0.2158 0.1184 0.181 0.2608 0.1634 0.0617 0.0621
## ar9 ma1 ma2 ma3 ma4 ma5 ma6 ma7
## 0.1999 -0.2992 -0.5366 1.1572 -0.3980 -0.3814 0.1207 -0.2981
## s.e. 0.0492 0.1703 0.1518 0.1798 0.2069 0.1867 0.1476 0.1286
##
## sigma^2 estimated as 348471: log likelihood = -5698.81, aic = 11431.62
#Affichage des résidus
tsdisplay(ts_cnt_season_adj_arima$residuals, lag=20)
On constate que si ce modèle optimisé a un AIC plus petit que le modèle auto-ARIMA (respectivement 11431.73 et 11496.75), la différence n’est pas grande et ce modèle n’explique pas mieux les résidus. En revanche ce modèle est beaucoup plus complexe que le modèle auto-ARIMA de paramètres 1,1,1. Il est donc préférable de continuer avec le modèle auto-ARIMA pour éviter l’overfitting.
#Prédictions sur 30 jours avec le modèle auto-ARIMA
plot(forecast(ts_cnt_auto_arima, h=30))
#Original
plot(ts_cnt_season_adj, col="red")
legend(x="bottomright", legend=c("Original", "Fitted"), col=c("red", "blue"), lty=1:2, cex=0.6)
#Forecasted
lines(fitted(ts_cnt_auto_arima), col="blue")
Les prédictions sont plutôt fidèles. Les gros pics ont été lissés.
end_time = time(ts_cnt_season_adj)[700]
train_set <- window(ts_cnt_season_adj, end=end_time)
test_set <- window(ts_cnt_season_adj, start=end_time)
#Initialisation du modèle ARIMA (9,1,7)
manual_fit <- Arima(train_set, order=c(9, 1, 7))
manual_forecast <- forecast(manual_fit, h=30)
#Calcul de la précision du modèle
print(paste("Précision (RMSE) du modèle manuel: ",
accuracy(manual_forecast, test_set)[2,"RMSE"]))
## [1] "Précision (RMSE) du modèle manuel: 602.038169309621"
#Initialisation d'un modèle avec auto-ARIMA
auto_fit <- auto.arima(train_set, seasonal=FALSE)
auto_forecast <- forecast(auto_fit, h=30)
#Calcul de la précision du modèle
print(paste("Précision (RMSE) du modèle auto-ARIMA: ",
accuracy(auto_forecast, test_set)[2,"RMSE"]))
## [1] "Précision (RMSE) du modèle auto-ARIMA: 694.107740668932"
#Affichage et comparaison des résultats
plot(ts_cnt_season_adj, col="black", main = "Comparaison des différents modèles")
legend(x="bottomright", legend=c("Données", "Modèle manuel", "Modèle auto"), col=c("black", "red", "green"), lty=1:2, cex=0.6)
lines(fitted(manual_fit), col="red")
lines(fitted(auto_fit), col="green")
D’après la mesure RMSE et comme l’indique le graphique, le modèle ARIMA de paramètres (9,1,7) minimise d’avantage les erreurs que le modèle manuel.
#Rappel de la précision du modèle manuel
print("Précision du modèle manuel:")
## [1] "Précision du modèle manuel:"
accuracy(manual_forecast)
## ME RMSE MAE MPE MAPE MASE
## Training set 28.95805 598.7046 435.867 -1.568717 10.99067 0.1815598
## ACF1
## Training set -0.002659762
#Rappel de la précision du modèle auto-ARIMA
print("Précision du modèle auto-ARIMA:")
## [1] "Précision du modèle auto-ARIMA:"
accuracy(auto_forecast)
## ME RMSE MAE MPE MAPE MASE
## Training set -2.446428 635.8652 455.2605 -2.656829 11.65282 0.1896381
## ACF1
## Training set 0.005113433
#Rappel de la précision du modèle Holt-Winters
print("Précision du modèle Holt-Winters:")
## [1] "Précision du modèle Holt-Winters:"
accuracy(ts_cnt_HW_arima)
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 2.77562 1009.124 740.1123 30088.77 30167.76 0.5646584 -0.009851498
La comparaison des scores nous encourage à utiliser le modèle manuel, plutôt que le modèle auto-ARIMA ou Holt-Winters. Par exemple le RMSE du modèle manuel est d’environ 599 contre 636 pour le modèle auto et 1009 pour le modèle Holt-Winters, ces derniers se trompent donc davantage.
#Prédictions sur 25 jours avec le modèle auto-ARIMA
plot(forecast(auto_forecast, h=25))
#Prédictions sur 25 jours avec le modèle manuel (9,1,7)
plot(forecast(manual_forecast, h=25))
#Test de Shapiro
shapiro.test(manual_forecast$residuals)
##
## Shapiro-Wilk normality test
##
## data: manual_forecast$residuals
## W = 0.95935, p-value = 5.265e-13
#Test de Ljung-Box
Box.test(manual_forecast$residuals, type="Ljung-Box")
##
## Box-Ljung test
##
## data: manual_forecast$residuals
## X-squared = 0.0049733, df = 1, p-value = 0.9438
Le résultat du test de Shapiro indique que les résidus ne suivent pas une distribution normale (la p value est plus petite que 0,05). Toutefois d’après le test de Ljung-Box, les résidus ne sont pas auto-corrélés (la p-value est plus grande que 0,05). Comme vu précédemment il y a encore de l’information contenue dans le bruit qui n’est pas expliquée par notre modèle.
Bien que le RMSE du modèle manuel soit meilleur que celui du modèle auto-ARIMA, le modèle manuel est en revanche beaucoup plus complexe puisqu’il est de paramètres 9,1,7 contre 1,1,1. L’amélioration du RMSE apporté par tous ces paramètres supplémentaires ne vaut sans doute pas le coup par rapport à la complexité ajoutée et aux ressources de calcul nécessaires. De plus ce modèle risque plus probablement l’overfitting. Il vaudrait donc mieux utiliser le modèle auto-ARIMA.
Toutefois avec ces deux modèles, une partie de l’information reste contenue dans les résidus. Afin d’améliorer notre modèle de prédiction nous pourrions tester un modèle permettant de prendre en compte plusieurs saisonnalités (ici nous avons au moins des saisonnalités de semaine en plus des années). La fonction tbats() du package forecast le permet.
Une autre option serait d’inclure d’autres variables. Ceci peut être fait avec une régression. Une régression dynamique permettrait également de prendre en compte l’erreur et de la considérer comme une variable à expliquer et non comme un bruit comme dans une régression classique. La technique de stepwise regression pourra ensuite être utilisée pour sélectionner le modèle offrant le meilleur équilibre entre complexité et précision.